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We show that all perfect odd integer squares not divisible by 3, can be 
usefully written as \fN = a + 18p, where the constant a is determined by 
the basic properties of N. The equation can be solved deterministically by an 
i— i efficient four step algorithm that is solely based on integer arithmetic. There 

is no required multiplication or division by multiple digit integers, nor does 
the algorithm need a seed value. It finds the integer p when N is a perfect 
OH square, and certifies N as a non-square when the algorithm terminates without 

a solution. The number of iterations scales approximately as log${^/N /2) for 
square roots. The paper also outlines how one of the methods discussed for 
squares can be extended to finding an arbitrary root of N. Finally, we present 



a rule that distinguishes products of twin primes from squares. 



1 Problem and Summary 

a 

In cryptography and broader applications, such as the general factorization of in- 
tegers, it is often of interest to determine quickly if an integer A is a true square 
or is a product of two closely spaced integers, such as twin primes. Numerous non- 
deterministic methods of varying efficiencies are currently applied to the problem 
of finding the square root of an integer. These methods include digit-by-digit cal- 
culations, Taylor series, Newton's method ref [1], Pell's equation ref [2], and even 
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sophisticated hardware implementations ref [3]. For very large integers, these meth- 
ods, may not yield rapid conclusive results. In contrast, the algorithm presented here 
provides a conclusive answer after a specific number of steps, determined by the size 
of the N. Throughout this presentation, all variables and computations are based 
on positive integer arithmetic alone. 

The number of required iterations to find the square root of N is approximately 
log 9 (y/N /2) — 3/2. The same number of iterations is required to certify that iV is a 
non-square. 

2 Properties of N 

Suppose the object is to find the square root of a general integer M. The presented 
algorithm works by focusing on N, the odd factor of M that is not divisible by 2 or 

3 to any power. 

M = 2 k *3 l *N (1) 

Clearly, if either k or I is odd, M is not a perfect square. If k and / are both even, 
then M is square if N is square. Since factors of 2 and 3 can quickly be identified 
in even large numbers, we can choose a representation for N which do not contain 
these divisors. 

One such representation is the reduced residue system [/ 18 with the elements 
{1,5,7,11,13,17} or {±1,±5,±7} modulo 18, where each element [r] is a residue 
class in Uig. By definition, ref [4], all [r] are relatively prime to 18 from which it 
follows that g.c.d. (r, 18) = 1 and the congruence r*N = 1 (mod 18) has a solution. 
It also happens to be true that all [r] in C/ 18 are relatively prime to each other, which 
is not generally true for reduced residue systems; but is true, for example, for systems 
modulo 3, 6, 8, 10, and 18. 

Since the integer 18 only contains the divisors 2 and 3, Uis contains no even 
numbers and no odd numbers divisible by 3. Therefore, U\g is an ideal representation 
of all possible iV-values because every iV defined in equation 1 can be uniquely 
assigned to one of the six elements {1, 5, 7, 11, 13, 17}. According to Euler's 
theorem, for Vr relatively prime to 18, r 6n = l(modl8). See figure 1. 

To test for divisibility by 3 one traditionally computes the sum-of-digits (SOD) 
of N repeatedly until a number between one and nine results, often called the Digit 
Root (DR). If DR is 3, 6, or 9 then N is divisible by three. This test is equivalent to 
computing N(mod9) where the latter calculation involves division by 9 as opposed 
to summation of digits. 
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Here, we extend this test to cover N(modl8) replacing division by 18 with the 
addition of the digits in N. We will show that a mapping exists between the elements 
in U 18 , N(modl8), and DR(N), N(mod9), of the form: 



[N{modl8)](mod9) = N(mod9) 



(2) 



Proof : 

For c < 18 write N(modl8) as N = C+18X = c+9(2X) and for d < 9 write N(mod9) 
as iV = d + 9Y. It follows that (c + 9(2X))(mod9) = c(mod9) and (d + 9Y)(mod9) = 
d{mod9) = d. Therefore, the equality in equation 2 is satisfied when: 



This condition is uniquely met for all N in Ui$; and therefore the DR method 
can be applied to the elements in Ui$ as long as it is recognized that we have to use 
c = 11 for d = 2, c = 13 for d = 4, and c = 17 for d = 8. For 1, 5, 7: c = d. To 
find which element an arbitrary N belongs to, we apply the condition in equation 
3. By doing so, we have substituted division by 18 with addition. For example, if 
N = 512346251, SOD(N) = 29 with DR(29) = 2, then N belongs to element 11, 
written as [11] or 11 + 18X. 

3 Determining solutions for a for known N 

Writing any integer in Ui$ as a + 18p, with a is uniquely determined by DR(N), any 
square in U 18 can be written as: 



From this expression DR(N) equals DR(a 2 ) because the remaining terms contain 
the factor 18, which does not affect the DR. Since a is limited to the values 1, 5, 
7, 11, 13, 17, the squares in U 1S are therefore found only in: [1], [7], [13]. This 
simplification comes with a computational cost for squares because two a-values are 
possible for each element: 



c(mod9) = d 



(3) 



N = (a+ 18p) 2 = a 2 + 36ap + 18V 



(4) 



• [1] for a = 1, 17 

• [7] for a = 5, 13 



• [13] for a = 7, 11 
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Figure 1: The complete modl8 power multiplication table for U\g. Example: a = 5 
to the 4 th power belongs to residue class [13] because 5 4 = 625 and N is therefore 
congruent to 13(modl8). Conversely, N can't be a power of 2 + 6n, and specifically 
a square, if Nmodl8 has a value other than 1, 7, or 13. 



Therefore, two parallel tests are needed for each square root calculation. 

As argued above for squares, the value of N{modl8) is required to determine the 
possible values of a for any root. From the first, third, and fifth column in figure 1, 
all even powers belong to the residue classes [1], [7], [13]. Also note that six tests for 
the constant a are required for powers of 6n, three tests for powers of form 3 + 6n, 
two tests for powers of 2 + 6n and 4 + 6n, and only one test each for powers of 1 + 6n 
and 5 + 6n. All such tests are independent and can be performed in parallel. 

4 Algorithm Description - Squares 

The algorithm is divided into three parts: 

1. Setup 

2. Initialization 

3. Loop 

4.1 Setup 

With and a given rewrite condition 4 as: 

N -ap 2 
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P 2 (5) 
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Here N = ^f 2 - is guaranteed integer when (DR(N),a) assumes any of the six 
combinations (1,1), (1, 17), (7, 5), (7, 13), (13, 7), (13, 11), see figure 1. Therefore N 
is square with a root of the form (a + 18p) if there exists an integer solution for p 
satisfying condition 5. We will iteratively test for the existence of such p and for this 
purpose define the unique representation of p as: 



p 



(6) 



i=0 



with 1 < b i+ \ < 9 and where i-max is the smallest i-value for which the termination 
criterion: 

No -apt 2 / 7 x 
9 ^ Pi ( ? ) 

is met. When the equality sign is valid, then N is square; and when the reverse is 
true a solution may still exist. In this case the iterations must continue until the 
termination criterion is met. 



Pi 



h + 9b 2 + 9% + ... + 9% +1 = + 9^+1 



(8) 



4.2 Initialization 

Set % = 0, bo = 0, and p-\ 



0. 



Initialization 


N 


frac 


b 


Po 


fo 


i = 


N-a A 
3fi 


N — a*bi 
9 







h *Po 



Figure 2: Init shows the algorithm initial conditions, where b\ is determined from 
the condition that fraco be integer. 



It is useful to define the quantity f = bi * bi = p * b\ and call the left hand side 
of equation 7 fraco. 



4.3 Loop 
4.3.1 i = 

Using equation 7 and the initial values from figure 2 condition 5 becomes an evalu- 
ation of: 

g == h * h (9) 



We know from the unique representation of p that there exists one and only one b\ 
value for which frac is integer and that b\ must satisfy DR(N ) = DR(a*bi). This 
condition can be implemented in a nine-entry look-up table. With b\ now known, 
condition 7 can be evaluated. The evaluation terminates when fraco < fo, continues 
for frac > f , and declares a solution for frac = f . 

4.3.2 i = 1 

In the next step p\ — b\ + 9b 2 . Inserting p 1 and pi into condition 5 yields: 

frac - f - ab 2 
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= b 2 (2b l + % 2 ) (10) 



Setting fraco — fo = JVi and realizing 2b\ + 9b 2 = Po + b\ + 9b 2 yields a condition 
similar to 9; but now with b 2 as the unknown quantity: 

g ==& 2 (P0+Pl) (11) 

As before, b 2 is found from the integer condition of the left hand side, frac\. As i 
increases, the number of terms on the right hand side, fi, increases as well. However, 
since each f\ can be expressed in terms of the previous and current value of pi, there 
exists a general recursion relation for f as shown in figure 3. 

4.3.3 General Case: i > 1 

One can express the i'th value, for i > 0, of iVj as: 

N _ = N l _ l -a*b l -V*b 2 l -2*b l * Efclp" 1 9 fc * b k (12) 

and 

Ni ~ a * bi+1 = 9* * b 2 +1 + 2 * b l+l * g 9*" 1 * 6 fe (13) 
y k=i 

Using the definition of pi, equation 8, the right hand side of equation 13 becomes: 

fi = h+i * (9* * b i+1 + 2 * pi_x) = b i+1 * (pi + pi_x) (14) 
The algorithm logic flow diagram is shown in figure 4. 
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Figure 3: Illustration of the general recursion relation for all the variables as a 
function of i. Reading the tables from left to right gives the order of computation. 
The quantity bi+i can assume the values 1 through 9 and is uniquely determined from 
the integer condition: Ni ~ a * bi + 1 _ 



5 Positive and Negative Certification 

Here we show that N is a perfect square if and only if the condition fraci = fi is 
satisfied after % steps. Conversely, N is not a square if fracj < f) without having 
satisfied the condition for squareness above for % < j. Of general concern could be 
that the algorithm would map different N values, squares and non-squares, to the 
same p. As we shall see, this is not the case. 

Combining the square condition: N = (a + 18p) 2 with equation 6 yields the test 
condition for squareness: 

N == (a + 18(6i + %2 + ... + 9%+i)) 2 (15) 

As discussed in section 4, equation 15 is equivalent to evaluating: 

iV - a(6i + 9b 2 + ... + 9% +1 ) == 9(6i + 9b 2 + ... + 9% +1 f (16) 

Section 4 also showed how each 1 < b^ < 9 can be uniquely determined for each 
step in i, with frac^ decreasing and fi increasing with % until condition fraci < f\ 
is satisfied or superseded. At step % all values of b are uniquely determined in the 
expansion of p\ and because the representation of p is a unique expansion of bi mod 
9, we see the validity of the perfect square condition stated above. If iV tests to be 
square, its root is also unique per equation 15. 



6 Generalization to Cube and Higher Roots 

For cubes, iVo is determined from N = (a + 18p) 3 as: 

N — (7 3 

iVo = = a 2 p + 18ap 2 + 6* 18p 3 (17) 

3 * 18 

Recall, from figure 1 that three parallel tests for a have to be conducted for this 
case. The first equation to be evaluated becomes: 
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Figure 4: Logic flow diagram for finding the square root of an odd integer N not 
divisible by 3. 
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N ° ~f bl = 2ab\ + 126? (18) 
9 

Although two divides both sides, we will maintain the same nomenclature as for 
squares. In this case frac = (N — a 2 bi)/9 and f = 2b\(a + 6&i). Proceeding as for 
squares Ni = frac — f becomes: 

JVi = ^o-° 2 *&i-9*/o (ig) 
9 

We return to equation 18 and insert the next value for p; pi = b± + 9b 2 yielding: 



N - a 2 (h + 9b 2 ) 



2a{h + 9b 2 y + 12(6i + 9b 2 f (20) 



9 

As in the quadratic case the termination criterion is met when: /raq < fa. 

Observe that the same equations and steps applies to cubes as for squares. The 
difference between the two methods lies in the existence of a recursion relation for 
squares. Perhaps recursion relations can be found for higher powers. If not, the 
above step-by-step approach can be applied to all powers. Although the formulae 
for a given power contain increasingly more terms as % increases, they only have to 
be derived once per root. 



7 Number of Iterations Required 

The smallest possible value of N of the form (a + 18p) * (a + ISp) is N = 361 for a = 1 
and p — 1. Smaller squares are possible for p = but are considered trivial in this 
context since the focus here lies in computing roots of very large integers. In order 
to find an approximate formula for the number of iterations, i, we recall that fraci 
decreases by a factor of 9 per step, whereas fa increases by a factor of 9 per step, 
and that the algorithm terminates when the two quantities are equal or fraci < f\. 
Therefore, the relative size of fraci and fi after % steps are: 

f rac o f (9V] 
= fo * 9 (21) 

From equation 5 and the conservative approximation of a = 0, we find fraco = 
N/ (9 2 *2 2 ). Setting fo — 9 conservatively in equation 20 yields the desired expression 
between N and the number of iterations i, ^ = 9 2i+3 . This can be rewritten as: 

i = log*(^f) ~ I (22) 
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Expression 21 represents the number of iterations required for each a value and 
therefore doubles when two tests have to be performed. Recall from equation 1 that, 
iV does not contain any factors of two and three. If iV is not a square, the same 
number of iterations are required. 

8 Numerical Examples 

8.1 Square Root Numerical Example 

As a point example consider \fN for N = 1429822969 = 37813 2 with DR(N) = 7. 
The allowed values for a in this case are 13 and 5. We will test for a = 13 first and 
see how the algorithm works when a solution is found, in figure 5, and then show 
what happens when there is no solution by testing for a = 5, figure 6. 



i 


a = 13 


DR(Ni) 


bi 


bi+i 


Pi 


fraci 


ft 




N = 1429822969 


7 















N = 39717300 


3 





3 


3 


4413029 


9 


1 


TVi = 4413020 


5 


3 


8 


75 


490324 


624 


2 


N 2 = 489700 


1 


8 


7 


642 


54401 


5019 


3 


N 3 = 49382 


8 


7 


2 


2100 


5484 


5484 


4 


iV 4 = o 




2 











Figure 5: Illustration on how the square root of N = 1429822969 can be found in 
four iterations. The algorithm terminates with a solution when fraci = /j. Thus 
p = 2100 = (3 + 9 * 8 + 9 2 * 7 + 9 3 * 2) with solution t/N = 13 + 18 * 2100 = 37813. 



8.2 Cube Root Numerical Example 

The following example shows the algorithm steps for finding the cube root of a 
thirteen digit number N = Q, 177, 847, 762, 549. Since DR(N) = 1, a can take on the 
values 1, 7, and 13. We will test iV for a = 7 since this is the solution in this example. 
The cubic root in this example is not computed by using a recursion relation but by 
step-by-step calculation of fraQ and f\. 
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Figure 6: Illustration on how a square root test of a ten- digit number N is concluded 
in four steps when N is not square. The algorithm terminates when fraci < fi which 
occurs for i = 3. 
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Figure 7: Example on how the cube root of a thirteen digit number N can be found 
in four iterations when the correct value for a is selected. The algorithm terminates 
with a solution when fraci = f. In this example p = 2 + 9*5 + 9 2 *3 + 9 3 *l = 1019; 
and the cube root of N is therefore: 7+ 18 * 1019 = 18349. 

9 Twin Primes 

Since U± s contains all odd primes, except 3, all twin primes are members of Ui$. In 
fact, since the difference between two twin primes is 2, twin primes can only occur 
between the residue classes [5], [7] (Type A); [11], [13] (Type B); and [17], [1] (Type 
C). The first 12 twin prime sets and their types are given in figure 8. 

Because the product between twin primes always belongs to [17] and squares 
always belong to [1], [7], [13]; the a, and therefore DR(N), in iV = a + 18p becomes a 
convenient differentiator between the two cases. 
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Figure 8: Illustration on how twin primes can be divided into three different types 
based on their residues classes (modl8). Note how the product of all three types can 
be written as 17 + 18 * X . 
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